{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Closed loop stability\n",
    "\n",
    "Stability of closed loop control systems appears to be easy to determine. We can just calculate the closed loop transfer function and invert the Laplace transform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy\n",
    "sympy.init_printing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = sympy.Symbol('s')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "K_c, t = sympy.symbols('K_c, t', positive=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are the systems from Example 10.4 in Seborg et al"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "G_c = K_c\n",
    "G_v = 1/(2*s + 1) \n",
    "G_p = G_d = 1/(5*s + 1)\n",
    "G_m = 1/(s + 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "K_m = sympy.limit(G_m, s, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "forward = G_c*G_v*G_p\n",
    "backward = G_m\n",
    "\n",
    "G_CL = K_m*forward/(1 + forward*backward)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALkAAAAlCAYAAAAeCUhkAAAABHNCSVQICAgIfAhkiAAABStJREFUeJzt22vIFFUcx/GPZhdvRBZZgfSYRXS/SBRlZkFUL7ILRdCFCLpAaUT3gsoXgYUFYWAEUtuLLlQSZUSWZm+yNMouRhmEGaRdzLKiwi724uyD88zOzM6zu+7O83i+MOzumTlz/r+dM2fOOXN+RCKRSAF74XtM6XUgKV7Ezb0OIlIt3sLTqbRr8AfuLMg3D0/tqKBymI7F2IBtuCjjmKOxGXt2Ma5IxfnZ9pZvdyzEJpxZkGcMfsGpHYqhhjkljjsH9+NC+ZUcPsANnQgsMvSZIlSW0zAJq/Ah+prku1hoLUek0kfgdqzFX/hB6D40o6ZcJU9SVMnvwzuDPN+wZGSvA6gAU4XKMkGo3J/jFHzdJN80obXclkq/DVfhehyKmXizc+GWZiVOwOgelF0pRvU6gAowFf/hBaH//VDJfH1CvzjN2XgNy+q/1+O99kJsiQ3YFQfgqx6UXxliSx4q+TJ8i+MHkW8PoTuS5iXchKVCn3ifnPx34/fEdllGWjv9/T/rnzt9Sx7hJ9wizEj8JrtfPAkvYzU+w2RhNub5nHNOEQayH2ELDss4ZgIOTmyLMD+V1qyCFvXJT6zvn9jkHJFhzmShIpxe/z0T/+DSxDG74dPEMXsKMzC3Yk2T84/Cr7ikRCw1nR14Xi08nXZ6dvY++dT65+r65ytCl+EJYeC5AhcIMy7L68dsqX8uwYNCd2RTPe0O4eXQKuFmuRJb8XYHYx4ntPL99OFYYabnm0T6dLzewXIjQ5S5sgdlTwpTf5OFOenrcvKvwKzE73vxhdAf3oRXcUzJWGrKteQzhBY8vdUSx4wWniAnlSw7spMzG4/Wv4/E3ol9Z+FL7NLtoJowC2/0OojI0GG8MCW4RujWTEvtvxEHdjuoJlwrzNFHIpFIJBKJRCKRSM/pX0GXXmQUiUQikUhrVNFOVmQlq2K8ZcnTNRw1DaBVG1inqKKdrMhKlo73LrwvvG38sX7eIzscbxnaschlXYM5Gt+wfte5cEvRsqb0UtvjBCMA221gc3EeHmghsJryi47GCDfUwhbKaafcsfjYwNfzST7BOlyeSs+KdwYW4GScIaxfWSqsOGyXms5pIltX0TVYi/0T21ElY2lGTTldrWoaQKs2sCJqyl+YLDtZN61k5LcQWVayPPtbknH4F+cm0qqiiUZdeZrmaL7isrK2v2RL3qoNrFNk2cmqbCXLs78lGS/8x5sTaVXRRKOuIk0HCUt31+FZjY1fVXQ1XKvkUttWbWCdok+jnazKVrI+2fa3JI8IxolkzFXRRKOuPtmaVuIKYTHaRGHs8S6OsP0GroquQtvfUmGN9Ho802IBafvW38J66jKWriV4PJU2W3jcD9ZKNphyk+Q9Bg+p70sOIrPiTTIPGw1c+011NNGoq5mmfsYKMzDJmYxu6hqMpgG0agNL0o6lK89O1i0rGfl/XpaVrMj+9rBQCQ7P2V8FTTTqKtKUZjkeS6VV2vbXjg2siJryg4pmdrIdbSUj/8/LspLlxTtfcQVP0ktNNOoqY+kjmLg3CiaRLCpl++vvk7djA+sUaTtZN6xklLOTZVnJsuxvC4Tpq/Pr+ferp/c/hqukiUZdWZoI47PF9bz74h6hy9I/nz4kbH/t2sDyqBncXZq0k3XDSkZzO1mRlSxtf8s6z7ZELFXRRL6utCZ4ThjQbRVayUUGPqmGje2vyAbWKapoJyuyklUx3rLk6RqOmkrTzAbWKapmJ2tmJatavGUp0jUcNUUikUgkEolEIpEe8T/qIQiYw/VO/wAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\frac{K_{c} \\left(s + 1\\right)}{K_{c} + \\left(s + 1\\right) \\left(2 s + 1\\right) \\left(5 s + 1\\right)}$$"
      ],
      "text/plain": [
       "           K_c⋅(s + 1)           \n",
       "─────────────────────────────────\n",
       "K_c + (s + 1)⋅(2⋅s + 1)⋅(5⋅s + 1)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.simplify(G_CL)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = G_CL/s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now for the moment of truth. Uncomment the next line if you have a lot of time on your hands..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#yt = sympy.inverse_laplace_transform(y, s, t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So that didn't really work as we expected. Can we at least calculate the roots of the denominator?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPkAAAAaCAYAAAB1n/q8AAAABHNCSVQICAgIfAhkiAAABpRJREFUeJzt22usHVUVwPEfRQERigJKE8EXIG0D+CIIYs29gijyaImK0fg4GsT4IkaNophYNCq+akFjsFFEwA8qig+soVBNRCwiIjYIWgSvQLRqqVpqhNL2+mHtyRnnnrn3njmPObfOP5nMOfsxs/aavfesvdYeGhoaGhq64u1Yjy3pWIdTa5VodPgAfin08nf8AEfWKlFDQwWW4mU4HM/Ax/AIjq5TqBHhWrxRDOyjcDU2Yv86hWpo6Aeb8Za6hRhB9sEOnF63ILsy8zqkXYdJjHfIOwg/Tvk/Mhoz8CvwedwgzMBJXDmLegfjUvwZD2MCK/H4Psq2O14tOvPP+3jd2dKtblqpzHTHjj7Kt6/og5v7eM1uOBVrcD/+g3vwLRxfkzzdULXfgwewE/ML6UvEgNiB5TpPEHVwm2jgg7jT7Bp7KP6ayn4XF2pPXr/DAT3KdBS2Yjv+Icz3OuhWN88Sz7bTsTbVv6aP8n0DvxaT4bD5pGjPJnxZ9IGrsE30/9fWIFM3VOn3iM4/id8X0t8j1pWb8JK+idmmle47VqHuuFj/7pbqz6ax16Zy7yykr0jpl1SQI88eOAzH4BNCb1UdTC3D1U0Z61L9MyrWL/Jp/EXoqQot1fWyQLysNuKJhbzxdN17Kso1W1qqy08Pz/ZVhcLz8Z2U9gs8uaJAM9HSW4Mzxszc2KenMn801RrZV7yB/43Hdqi7ROjjbjyEv+FmfHwGua4Xb4sqtAxPN2Ucmerer/yt241uPissqcUVZMloqa6X56W63yvJ3yLekJ2o2geKtPTnuTLDsy128mPS+RbhDb4FZ4o32xLc2weB6uZF6bxGmGV5HsSN2BvHFfI+iJ/iucJ0XSFCQHvipTPccx72qi5y7WROw6/ovCbvRjcXC1N4HHcMQthZcJcwy4/FgYW8F4rJ/voO9XrpA7XxqML/bJAvFjPTJF6nuok3ihyRzhtK8u/CySL8tTalHYSP4Gc4UXSQPPmOciF+iPtEZ3mNmGnnaqz8MWJQ7tTZGulGN19M11omnG0LUvrWdAyLzXi/GKR3CL/MA2K5eoZwPhejId20c6TID/Ld8Oz0+83CZD0Otw9bqAGzXzr/qyQ/S39cLm2hMFM3mPpwiTV3xgIxKS5I11qPU4QfYC5yltBFNnEV6UY3b03ntYUyFwjn3jBZKSIql4r+nvEHXCbM8DzdtHNkOULb6far9PtDA7jPhJnDNPnjsi6uPWbmdeeqVObskvzMgjkvl3Yg/pnSvy98F/0MtWVMqFc3nbgx1SuLZQ9DNxP6r5f3iejHCuGn2RvP0XbKfqpQvpd2DkL+PGOmebb5N3lmqq8T2w9vxkeFg+rrXd50Olb637ckEbpZiq8JheS5rY/3pv2m3q8kf36hHDFLvwAfFuGw08Xa9Dqcj1v7JFvduimyGM8XDrfVJWWGoZt+62VMhNCuxrtz6bcKH9QGEVG6RNvL3ks7R+a5fk7MBu9I/58pHFEPC2fEIGkZngf57FTmSyX52Ux+Ykn+HjgJ39SOs+5ZQdbZ0lKfd/2iVGf5LMsPUzct1fXyGZ1DqBlZROnlJfn9aGdLDd71vGcdfiNMkt2FY2JhyQ0OFubFRhFSuB0v7kHgQfOTdD5Z5xDaCWL3000l9bcJz+tZwglzgHDK7GrsJZyuO4VXfTbMFd1kA/IJJflZeqe1d5Y+F9qJdiefJ8yH7WJwZ6zGu8TaY7WpGwcOEYNhf7HN7mixximLMY4Cd4vw2VPFF2N5LhDx8cuF45FwRh7a4TqHifjxvcKc3dV4pfZz7+RwY+7q5oZ0PgdPKuSdIib6h7S3Is/VdqK9Jl8k9levF2+xPF8QO2vOFQ6H8VyZVfiTCIlkMeey0NSgWJYO2iGZ47WdF5vw3kKdt4kHeLEwy+8UGyTGhfzn58qeizcIH8Vvhdf1ado7v95karx9VKiim4xz0nnVNNefq7q5SryJTxLPPvsabhFOE5Gm80RYjdFsZ9fP9vXCpi8zy+aJAT6Jb6f/T0n/j+2DwC3V1yfLTe+pnCipdwi+KrZWbhOT1UWmfnSzDFeIqMOWVHZCxIwPryBvt7QMXzeLUv59pt9XXqduWnpb0z5aWKk3Cdm3i8F7jVjK5RlEO1t6k3+5as+2K5YKxYzKhyoNDQ0d6GWAPpLq79MnWRoaGgZAL5/4bRbrtoXan2eeKUyZ4m6hhoaGOcoJYkfUVvHd9BpTPfANDQ0NDQ0NDQ0NDQ0NDQ0NDQ3/d/wXuNJDzrDAA18AAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$K_{c} + 10 s^{3} + 17 s^{2} + 8 s + 1$$"
      ],
      "text/plain": [
       "          3       2          \n",
       "K_c + 10⋅s  + 17⋅s  + 8⋅s + 1"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ce = sympy.denom(G_CL.simplify())\n",
    "ce.expand()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "roots = sympy.solve(ce, s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, that approach works, but is limited in the analytic case to 4th order polynomials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWsAAAD5CAYAAADhnxSEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X98k/W9///HlaRJm6Zpkv5Mf0ChLVAoUKGKzo0hgqhjMKcTnJu4qXDUOeePM9k8euOzIw533L7uzM/RMXUyzyYTPgJOURAUdfij4ijKD7EIhf5I27RN+itpkzTX94+UCLaFIm2TlNf9duuN9rqu5Ho1hqdvXnlf70tRVRUhhBDRTRPpAoQQQpyehLUQQsQACWshhIgBEtZCCBEDJKyFECIGSFgLIUQMkLAWQogYIGEthBAxQMJaCCFigC7SBYjIUhSlFPgGkAV4gb3ANlVVmyNamBDiJDKyPkcpinKjoij/An4BJAAHgQbg68DriqKsURRlVCRrFEJ8QUbW565E4GJVVb197VQUpQQoBI4Na1VCiD4pspCTEEJEP2mDnON62h2WE362KoryTCRrEkL0JmEtpqiq6j7+g6qqLuC8CNYjhOiDhLXQKIpiPf6Doig25LMMIaKO/KUUvwXeVRRlPaAC1wIrI1uSEOLL5ANGgaIoE4HZgAJsV1V1/wn7rD2tESFEBElYi1NSFOVfqqpOi3QdQpzrpGctTkeJdAFCiDPvWcsw/Bxz3nnngfx3F+JsnfWgR0bWQggRAySsRS8+ny/8vXymIUR0kLA+R/3nf/5nn9tbWlq47LLLwj9v3759uEoSQpyChPU56p133uH+++8/aVtdXR0zZ85k9uzZ4W02m224SxNC9EHC+hz10ksvsWfPHu6++24AKioq+PrXv85tt93Ggw8+GOHqhBBfdqbzrKWBOYL4/X4WL15MXFwc7733Ho899hhXXXVVpMsSYiQ669kgEtbnqN/97ndAKLB/85vf8I1vfIOZM2eG9x8fcQshBsVZh7WsDXKOamtrC3//05/+tNc2IcTgCAaDaDRn33GWkbUQQvSj0xfA0dpJQ2sX9W2dOFu7aOrw0er14/b6ae7w0dbpp6Orm6Cq0ukP0trpx98dpDuoElTh4vwU/nrLhTKyFl/NQw89xO23347Vau1z/xtvvIHH42H+/PnDXNlXc/nll/Paa69FugwRRYLBIK4OP47WTupaOmlo68LZ3omrw0eLx0+Tx0ebN0B7VwCPrxuNAs0dPvzdKoFgkGDP0DROq+DvPvU4VVEgIU6LXqtBr9VgjteRoNeSqNdxcWHqoPw+EtbnqMmTJzN//nzi4+OZNm0aaWlpdHZ2UlFRQXl5OXPmzOGXv/xlpMscsMbGxkiXIAZZMBiktTNAjduLw+2lsa2LurYumtp9uDw+Wrx+Wr2hUWyLN0CHL0CXPxge1Q60DaAAWo1CUrwOnVaDKV6DUa/FqNeRFK8j3WQgXq/FatSTatKTaoonLUlPhjkee3ICpvjhiVEJ63PU+vXr2blzJ7/5zW9IT0/H4XBgNpv5wQ9+wOrVq0lISIh0iWIECASCOFo6Q4Hb0klje6id0NTuw90TuL5AkGaPH68vQFcgSDCoEuj5OpHJoKW9q7vXOfRaBY1GQasoJOi12OL0GPVaTAYtqSYDRoMOq1GPLVFPWpKBNJOeTHM82VYjyQm6QeknDwcJ63PURx99xNGjR/nrX//Km2++edI+r9crYS3CAoEgtS1eat2dOFo6qWv10tTuw9nehdvjw+314w8Eae7w4fV30xUIggr+oEr3lwJXp1F6hbBBp6E7qBKn1WCI05Co12FO0GHQaTEn6LAY9diMejLNemyJ8aSbDWRZErAnx5MUHzecL0VESVifo/7t3/6Nyy+/nMOHD1NaWhrerqoqiqJw+PDhCFYnhkJ7Z4Bjrg5qXZ3UtXipa+3C5emiobULt9dPizf0QVl3MEiL14+/Ww23E/ob1R4Xpw19fqbXhcLWmqjHZNCi12mxJMRhMcaRajKQYdaTbk4IjWwtCaSa9DEzso00mQ1yjrv11lt54oknhuS5X3vtNe688066u7u5+eabWb58+Un7u7q6uOGGG/joo49ISUnh73//O3l5eQD8+te/5umnn0ar1fLf//3fzJs375TnKi0tZdeuXUPye0Sbjq4Ax5o91Lg81LaEPjxztnXR3hWgvrUz1L/t8uP1dwMK7Z1+TvP5GPBF79Zk0KHRKCTEhVoJ5oQ47MnxxMeF2gppSYaeNkICOckJWE16FEWWPT8NmQ0izs5QBXV3dze33347r7/+Ojk5OZx//vksWLCAiRMnho95+umnsVqtHDp0iLVr13Lffffx97//nf3797N27Vr27dtHbW0tc+bM4bPPPkOr1Q5JrZESDAZpaOviaJOHapeHutZOHO5Omjp8NPXMWPB1d+Py+On0d+PvVgme4oMzg05DVyCIooTaDQadBkuCDqsxDpNBh8UYR1qSgeQEPek9gWtPjicnxYg9KR6dTka40UzCWgyJsrIyCgoKGDt2LACLFy9m06ZNJ4X1pk2bWLFiBQDXXHMNP/nJT1BVlU2bNrF48WIMBgNjxoyhoKCAsrIyLrrookj8KgMSDAZxtvuobGynqtlLTYuX+pZOnO1dtHj8NHv8tHX68fiO93RVfKebDkYogDUaBaNeR2LPjARTvI60JANWo550swG7OYFsazw5lgQyk+OlrTBCnVEb5PLLL1eHY4qU0+kkLS1tyM8zWKTe3lwuF62trYwePRqApqYmOjo6GDVqVPiYffv2UVhYiF6vB+CTTz6hqKiI2tpaTCZTeMW/Tz/9lIyMjF5zwp1OZ3jKXiAQYPLkyYNWf1BV8QWCoa/u0Fd3txru43arKkE1NNKFgfUHFUBRFDQK6DQaNApotQo6jYY4rUKcVkOcVoNeF5qrO1ydBXn/Dr2PPvpon6qqxWfzHFHZs461/qPU29u6devYsmULTz31FADPPfccZWVl/OEPfwgfM2nSJLZs2UJOTg4A+fn5lJWV8eCDD3LRRRfxgx/8AIDU1FT++Mc/cvXVV/d7vtP9Th5fgMPOdiobPVS5PNS6QxdHOFq9uD1+2jpD83QD3aEpY1+exfBlOk1omphBpyXVpMcQpwnNWEiOZ+Pa53jgnjsYZUsg15pAVnICcXHR28KR9+/QUxTlI1VVS09/ZP+kDSKGRE5ODlVVVeGfq6urycrK6vOYnJwcAoEALS0t2Gy2Xo/1+/29HhsMBql2eTnU0MGx5g4cLZ385G//wu3xUd/aRae/G7fHT2egO3zZ75dpNdAdDF19FqfVEK/TYDPrsSboMfbM0U03G7AnxzPKmkiuLYG81ESM+lP/tXn94Tf5/oz/+gqvmhD9k7AWQ+L888+noqKCI0eOkJ2dzdq1a/nb3/520jELFixgzZo1XHTRRaxfv55Zs2fzeUM7qVMv5ber/4LnxT0cbnCT+J0V3PeOl/ZtW/F3q3R0BXqFb2N7Fy9/7ABC7QajXotWo5BqMpCcEIfVqMeeHE+qSU+21UiuzciYlARyrYnywZqICVEZ1kuXLo10CWdE6u1Np9Px+OOPM3fePAJxJi6/9iY+cBm55xdPYbKl052YQlNHEbWmTPLu3QAaI4xdwpz/7+3QE0y4gufKqgHQZxZS7fJi0GmxGONISdSHwtcST0qigWxLAr/9h5GX7vw6Y1ISiT/NyHeoyfthaMVavT1Wn+0TRGXPWkS/QCDI4aYOPq1r7Znz66Wq2YOzzYfXH6Cpw0cwSM9c35Ml6rV09CycY9BpMeq1pCTqQ6NgYxx2Szy5ViNjUo3kpyWRlWw47bS9WOxjinOKzLMWgysYDHKkycNBRxsVzjaqmkMfxjW2d+Hy+PEHumnt7N2GOJFRryVOq5BqNmDUa7Enx5OeFB9qPaQmUpCeSF5KokwxE+IMSFifQ9o6/eyrbeGz+nY+d3ZQ3Ry6EKO5w4cvEMTt9fc7C0IhdCmx1ajH2tOGyEwOXTI8JjWRgjQT4zKTsBj1w/tLCXGOiFhYV1VVccMNN1BXV4dGo2Hp0qXceeedJx2zY8cOFi5cyJgxYwD47ne/G9Gbuebl5ZGUlIRWq0Wn0/X6Z7eqqtx5551s3rwZo9HIs88+y7Rp04alNndHF3trW/m0ro3DznY+rW6kfO9+1MRUVJ0BVaNDUXqPZFVVRVGDKL4OTAl6xuZmkHl8FJyWSEG6kaLMZMwJgx/CP/7xj3n55ZdJT09n7969ADQ3N7No0SIqKyvJy8vjhRde6HPN7TVr1vDQQw8B8B//8R+DXttA6/33f/93/vGPf6DX68nPz+fPf/4zFoul12NP994ZrnpXrFjBn/70p/A85Ycffpgrr7yy12NPt1TAcNW7aNEiDh48CIDb7cZisVBeXt7rsZF4ffvLsIG+hxVFWQIcf/M+pKrqmlOdL2I9a4fDgcPhYNq0abS1tTF9+nQ2btx40hVuO3bs4NFHH+Xll18erNOelby8PHbt2kVqat+LiW/evJk//OEPbN68mQ8++IA777yTDz744KzP6/N1c6ChjX01rXxW30Z9aydHGjto6rlLRZc/2Od/mHidBhVINOhorKmkdFIhBXYrY1MTGZtmouXofv78P79jc4Re37fffhuTycQNN9wQ/sv585//HJvNxvLly1m1ahUul4tHHnnkpMc1NzeHe9SKojB9+nSSk5PZvXv3sNe7detWZs+ejU6n47777gPoVS+c/r0zXPWuWLECk8nEvffe2+/juru7GTdu3ElLBTz//PMn/d0crnpPdM8995CcnNzngC0Sr29/Gfbss8/29R4+qWetKIoN2AWUEsrVj4Dpqqq6+jtfxEbWdrsdu90OQFJSEkVFRdTU1Az5G2Iobdq0iRtuuAFFUbjwwgtxu904HI7w79mfFo+Pj2ta2FvTQkV9O3WtnVQ3e+jwdffZmjh+5wq9ToNJryPboifXmkC2NYHRKUbGZyRRnJ1MiskAhALl//yfZ3nxyZ0nPc8O56dEsms8c+ZMKisrT9q2adMmduzYAcCSJUuYNWtWr/DbsmULc+fODV/hOHfuXF5//fWI1HvZZZeFv7/wwgtZv379kNcxUH3VOxADWSpgKJyqXlVVeeGFF3jjjTeGtIYz0V+GDeQ9DMwDXldVtRlAUZTXgcuB5/s7X1T0rCsrK9m9ezczZszote+9995j6tSpZGVl8eijjzJp0qQIVBiiKAqXXXYZiqKwbNmyXlOIampqyM3NDf+ck5NDTU0NZlsa5VVu9ta4+bSunWPNHdS1hNYC9ner+LqDvc6lVRQsRh325HjSTAZyrUbGpidSZDdTnGUm22occN1r167luuuu63NfNL2+APX19eG/AHa7nYaGhl7H9PU6+/3+YauxP8888wyLFi3qc9/p3jvD6fHHH+cvf/kLpaWl/Pa3v+31T/S+Xt/B+Bfi2XjnnXfIyMigsLCwz/2Rfn1PzLCBvIeBbKDqhJ+re7b1K+Jh3d7eztVXX81jjz2G2Ww+ad+0adM4evQoJpOJzZs3853vfIeKiooIVQo7d+4kKyuLhoYG5s6dy4QJE5g5cybBYJDDjR00mcex5hMPTxx4nxq3l2Pn3cbi/1dPYP2WXs+l1Sgk6rVkW79Y27cwI4nJOWYm2c0YDYOzqLrP5+Oll17i17/+da990fb6DtQZtu6GxcqVK9HpdFx//fV97u/vvTPcbr31Vh544AEUReGBBx7gnnvu4ZlnnjnpmL5e30gvgfr888/3O+CAyL6+p8qwU+jrBT3lGzuiYe33+7n66qu5/vrr+e53v9tr/4m/+JVXXsltt91GY2PjsPaljgsEgjiDRl795xH217YQd+UvuWtLAx1bX6Uz0DMyHnUpjqN+oAmtAkEVci0GMpITyUs1MiEziSk5FoqzzMN24carr77KtGnTyMjI6LUvml7f4zIyMsKtI4fDQXp6eq9jcnJywv/MhNCl7HFxkbtjyJo1a3j55ZfZvn17v6F2/HL59PR0rrrqKsrKyiIS1ie+D2655ZY+b4g8kKUChlMgEODFF1/ko48+6veYSL2+fWXYQN7DhEbSs074OQfYcapzRSysVVXlpptuoqioiLvvvrvPY+rq6sjIyEBRFMrKyggGg6SkpAxZTcFgkIqGDsqONPFxdQsVDe00t/uoa+sk0B08eW6xmkiiVsFo0DEm1UCuzYimrZ69b/2DTU//nop9u/npT3/KW2VlQ1bvQJxqRDLcr+9AHL8Effny5axZs4aFCxf2OmbevHn88pe/xOUKfRazdetWkpOTh7tUIDRr4pFHHuGtt97CaOy7NdXR0UEwGCQpKYmOjg62bt0asVlNJ36GsmHDBoqLey8EN5ClAobTtm3bmDBhQnjBry+L1OvbX4YN5D0MbAEeVhTleA/qMuAXpz3hGXwNmnfeeUcF1MmTJ6tTp05Vp06dqr7yyivqE088oT7xxBOqqqrqH/7wB3XixInqlClT1BkzZqg7d+4clHO3ef3q63sd6sqX96k/+vMH6jce2a4WPfCqOmb5y+ro+07+mvAfr6rT/3OreuXvtqn5P3xIHb/wdnXcxVeqK371kKqq6kn1BoNB9bbbblPHjh2rFhcXqx9++OGg1PtVdXR0qDabTXW73eFtw/H6DtTixYvVzMxMVafTqdnZ2epTTz2lNjY2qrNnz1YLCgrU2bNnq01NTaqqquqHH36o3nTTTeHHPv3002p+fr6an5+vPvPMM+r06dMjUm9+fr6ak5MTfg8vW7ZMVVVVrampUa+44gpVVVX1888/V6dMmaJOmTJFnThxovrQQw8Nea391fuDH/xALS4uVidPnqx++9vfVmtra3vVq6qq+sorr6iFhYXq2LFjI1qvqqrqkiVLwu/Z46Lh9e0vw/p5D0No5sdTak+eAj8GDvV8/Ug9Tf6O6MvNm9q7eKfCSVllM586QlfjuXvuLXcijQJJ8XFkJRsYlZLIlJxkpo+yUZKbHPF1JsTAyOXmIsrJ5eYA7Z1+/nmokQ8ON/FxTQtHm7x0+AJ4fSevS6HXKmQkxZOZHE9RZhLT82x8LT+FdHN8hCoXQoiBibmwPljXxo6DDbx/uImKhnacbV34AidfFKLXarBb4slKTmByTjIX5Nn4WkHKadchFkKIaBW16RUMBik74mL7p/XhD/vcHh9BFeI0Cv6gikYBi1FPcXYiE+1mvlZg4+sFaSTFR25mgBBCDIWoCGtVVdlT5ebVvXXsrnLxWX07LR7/SaNlg05DjjWB8ZlmZuTZmDspg9EpiRGrWQghhlNEwrqpvYuX9tTy7qEm9lS7aWzvCk+LUwBDnIYxaYlMzU3m4rGpzJmYjsVoiESpQggRFYYlrPfXtrCpvJa3K5xUNnpOWpBep1HIsiQwLdfCxYWpXDYxA2uiBHMsG8iqY+Xl5dx66620trai1Wq5//77w5dq33jjjbz11lvhudPPPvssJSUlw/57CBFNhmTq3sfVbtZ9VMXOiiZq3d4vrvADzPE6JmWbuSDPxoKpWeSnJ51hySLaDWTlvM8++wxFUSgsLKS2tpbp06dz4MABLBYLN954I/Pnz+eaa64Z8Dll6p6IJJPJRHt7OxBaffPOO+9k+/btjBo16vgh0TF1r9bt5bn3jvJ2hZNDDW10BUKZrgApJj0X51iYNymDb02xkzhIa16I6DWQVcfGjRsX/j4rK4v09HScTmefa0ELESu2b9/OHXfcwdatW08M6kHxlcI6GAzyyicONpXX8v7hZtq7AkDo4hKrUc/McRa+PSWLyydloo879b3zxMgzwFXHwsrKyvD5fOTn54e33X///fzqV7/i0ksvZdWqVRgM0hoT0e2dd97hlltuYfPmzSe9l79MUZQM4ElgbM+mW1VVffd0z39GbZD/efOQumlPLYfq2+hWQzc+DapQZE/iW1PsLDo/F5OMnM8Jc+bMoa6urtf2lStXsmTJEtxud3ib1WoNr+PxZQ6Hg1mzZrFmzRouvPDC8LbMzEx8Ph9Lly4lPz+/z7UeVq9ezerVoZtGO51Ojh49Ohi/mhBnLC4ujqSkJHbs2MGUKVP6OiTcBlEU5e/Ae6qqPqYoihYwqaracrpznFFY5y1/RQVISzIwa1waP754NEVZ8s9WcbLx48ezY8eO8Kpjs2bNCt+a6UStra3MmjWLX/ziF3zve9/r87kGercg6VmLSDIajcyePZv8/Hx+//vf93XIiWHtBHJUVe06k3Oc0Y1CfnJJAeUPzuXD++fwX9+bKkEt+nR81TGg31XHfD4fV111FTfccEOvoHY4HEBo/v3GjRv7XBlOiGii0Wh44YUX+PDDD3n44YeH5iSnW+lJHaJV98TINZCV85577jlVp9OFVyubOnWqunv3blVVVfWSSy5Ri4uL1UmTJqnXX3+92tbWdtpzDseqe0L0JzExUVVVVW1qalInTpwYXjHwBOEcBdYCP+v5XguY1QHk74hedU+cO6QNIiLpxKl7VVVVzJw5k8cee+zEf1We2AbJAFYT+oCxm9AHjO+d7hwS1mJEkLAWUe6s51lH8ubWQgghBkjCWgghhsDKlSspKSmhpKQERVHKe77u/6rPJ20QMSJIG0REOWmDCCHEuUDCWgghYoCEtRBCxAAJayGEiAES1kIIEQMkrIUQIgZIWAshRAyQsBZCiBggYS2EEDFAwloIIWKAhLUQQsQACWshhIgBEtZi0DU3NzN37lwKCwuZO3duvzfL1Wq14VXJFixYEN5+5MgRZsyYQWFhIYsWLcLn8w1X6UJELQlrMehWrVrFpZdeSkVFBZdeeimrVq3q87iEhATKy8spLy/npZdeCm+/7777uOuuu6ioqMBqtfL0008PV+lCRC0JazHoNm3axJIlSwBYsmQJGzduHPBjVVXljTfe4JprrvlKjxdipJKwFoOuvr4eu90OgN1up6Ghoc/jOjs7KS0t5cILLwwHclNTExaLBZ1OB0BOTg41NTXDU7gQUUwX6QJEbJozZw51dXW9tq9cuXLAz3Hs2DGysrI4fPgws2fPZvLkyZjN5l7HKUrf67avXr2a1atXA+B0Ogd8XiFikYS1+Eq2bdvW776MjAwcDgd2ux2Hw0F6enqfx2VlZQEwduxYZs2axe7du7n66qtxu90EAgF0Oh3V1dXh475s6dKlLF26FAjdKUaIkUzaIGLQLViwgDVr1gCwZs0aFi5c2OsYl8tFV1cXAI2NjezcuZOJEyeiKAqXXHIJ69evP+XjhTjXyD0YxaBramri2muv5dixY4waNYp169Zhs9nYtWsXTz75JE899RTvvvsuy5YtQ6PREAwG+dnPfsZNN90EwOHDh1m8eDHNzc2cd955/O///i8Gg+GU55R7MIood9b3YJSwFiOChLWIcnLDXCGEOBdIWAshRAyQsBZCiBggYS2EEDFAwloIIWKAhLUQQsQACWshhIgBEtZCCBEDJKyFECIGSFgLIUQMkLAWQogYIGEthBAxQMJaCCFigIS1EELEAAlrIYSIARLWQggRAySshRAiBkhYi0HX3NzM3LlzKSwsZO7cubhcrl7HvPnmm5SUlIS/4uPj2bhxIwA33ngjY8aMCe8rLy8f7l9BiKgjt/USg+7nP/85NpuN5cuXs2rVKlwuF4888ki/xzc3N1NQUEB1dTVGo5Ebb7yR+fPnc8011wz4nHJbLxHl5LZeIvps2rSJJUuWALBkyZLwiLk/69ev54orrsBoNA5HeULEJAlrMejq6+ux2+0A2O12GhoaTnn82rVrue66607adv/99zNlyhTuuusuurq6+nzc6tWrKS0tpbS0FKfTOTjFCxGlpA0ivpI5c+ZQV1fXa/vKlStZsmQJbrc7vM1qtfbZtwZwOBxMmTKF2tpa4uLiwtsyMzPx+XwsXbqU/Px8HnzwwVPWI20QEeXOug2iG4wqxLln27Zt/e7LyMjA4XBgt9txOBykp6f3e+wLL7zAVVddFQ5qIDwqNxgM/OhHP+LRRx8dvMKFiFHSBhGDbsGCBaxZswaANWvWsHDhwn6Pff7553u1QBwOBwCqqrJx40aKi4uHrlghYoS0QcSga2pq4tprr+XYsWOMGjWKdevWYbPZ2LVrF08++SRPPfUUAJWVlVx88cVUVVWh0Xwxbpg9ezZOpxNVVSkpKeHJJ5/EZDKd8pzSBhFR7qzbIBLWYkSQsBZRTqbuCSHEuUDCWgghYoCEtRBCxAAJayGEiAES1kIIEQMkrIUQIgZIWAshRAyQsBZCiBggYS2EEDFAwloIIWKAhLUQQsQACWshhIgBEtZCCBEDJKyFECIGSFgLIUQMkLAWQogYIGEtBt26deuYNGkSGo3mlDcEeO211xg/fjwFBQWsWrUqvP3IkSPMmDGDwsJCFi1ahM/nG46yhYhqEtZi0BUXF/Piiy8yc+bMfo/p7u7m9ttv59VXX2X//v08//zz7N+/H4D77ruPu+66i4qKCqxWK08//fRwlS5E1JKwFoOuqKiI8ePHn/KYsrIyCgoKGDt2LHq9nsWLF7Np0yZUVeWNN97gmmuuAWDJkiVs3LhxOMoWIqpJWIuIqKmpITc3N/xzTk4ONTU1NDU1YbFY0Ol0J23vy+rVqyktLaW0tBSn0zksdQsRKbpIFyBi05w5c6irq+u1feXKlSxcuPC0j+/rRs2KovS7vS9Lly5l6dKlQOiGuUKMZBLW4ivZtm3bWT0+JyeHqqqq8M/V1dVkZWWRmpqK2+0mEAig0+nC24U410kbRETE+eefT0VFBUeOHMHn87F27VoWLFiAoihccsklrF+/HoA1a9YMaKQuxEgnYS0G3YYNG8jJyeG9997jW9/6FvPmzQOgtraWK6+8EgCdTsfjjz/OvHnzKCoq4tprr2XSpEkAPPLII/zud7+joKCApqYmbrrppoj9LkJEC6WvHuEpnNHBQgyX0tLSU87pFiLC+v7g5QzIyFoIIWKAhLUQQsQACWshhIgBEtZCCBEDJKyFECIGSFgLIUQMkLAWQogYIGEthBAxQMJaCCFiwBmFtatD7tghhBCRcEZhPevRHfztg2MEg3LVuRBCDKczCusJmUn8csMnXPXEu+ytaRmqmoQQQnzJGYX12qUX8ti0osagAAAT4ElEQVSiEmpcXhY8/k9+9Y99tHj8Q1WbEEKIHmcU1oqi8J3zstl+zzf54YWjOeBoY+Z/vcmf3j5Mp797qGoUQohz3lktkXrA0cojr33KjoNOsi0J3DtvHAunZqPRnPVqgEKcEVkiVUS5sw7FQVnPeuehRn796gH21rQy0W7mF1dM4OuFqf3eO0+IwSZhLaJcdKxnfXFBKi/d/nV+v7gEc4KOG/5cxnefeJe3PnP2eQNUMXKtW7eOSZMmodFo+g3PqqoqLrnkEoqKipg0aRK///3vw/tWrFhBdnY2JSUllJSUsHnz5uEqXYioNuh3ivEFgqz/qJr/++YhatxeSnIt/GxOId8clyYj7XPAgQMH0Gg0LFu2jEcffbTPu447HA4cDgfTpk2jra2N6dOns3HjRiZOnMiKFSswmUzce++9Z3ReGVmLKBcdI+sT6XUavj9jFG/eO4tff3cyzrYu/u+bh/jWf/+TTeU1BLqDg31KEUWKiooYP378KY+x2+1MmzYNgKSkJIqKiqipqRmO8oSIWUN2ublep+G6C0Khvag0l65AN3euLWfWozt4ducRPL7AUJ1axJDKykp2797NjBkzwtsef/xxpkyZwo9//GNcLlcEqxMiegz52iB6nYZrSnN5/a5v8qcbSskwx7PiH/u59X//xSOvfUqt2zvUJYhBNmfOHIqLi3t9bdq06Yyep729nauvvprHHnsMs9kMwK233srnn39OeXk5drude+65p9/Hr169mtLSUkpLS3E6nWf1OwkR7SJyd/Ndlc387YNjbCyvQVEULi/O5Edfy2P6aKv0tUeIWbNm9duzBvD7/cyfP5958+Zx991393lMZWUl8+fPZ+/evac9n/SsRZQ762DTDUYVZ6o0z0Zpno275o7jufePsrbsGNv31zMpO5mrp+WwsCSLRENEShPDQFVVbrrpJoqKinoFtcPhwG63A7BhwwaKi4sjUaIQUSciI+sv8/gCvPpJHX965zCf1rVhMuj4znlZXD9jNEV281CcUgyRDRs2cMcdd+B0OrFYLJSUlLBlyxZqa2u5+eab2bx5M//85z/5xje+weTJk9FoQp24hx9+mCuvvJIf/vCHlJeXoygKeXl5/PGPfwyH96nIyFpEuei4KGawqKrK7io3f33/GC9/XIstMY4McwLfK83h21OzMMfHDeXpRQyTsBZRbmSF9YncHh+v7a3jmZ1H+Ky+nfg4DVcU2/leaQ4XjrGFR2RCgIS1iHojN6zDJ1RVPq5u4YVdVbxUXku8Xku8TsOCkiyuOi+HgnTTcJckopCEtYhyIz+sT+T1dbPjYAPPf1jFPyucBFWYnJ3MVedl8+2pWaQlGSJZnoggCWsR5c6tsD5RQ2snL+2pZcPuGvbVtlI62kqCXsu3p2Qxb1ImyUbpb59LJKxFlDt3w/pEn9W3sf1APc+XVXGs2UOcVmFmYRrfnprF7AnpmBMkuEc6CWsR5SSsT3S8v/2PPbW88okDjaLgbO9iZmEqlxfbmVuUISPuEUrCWkQ5Cev+BIMq5VUu/vGxgy1766ht6USnUbgoP4Uriu3MnZhOWlJ8pMsUg0TCWkQ5CeuBUFWVPdUtvLrXwWt766hxeUiI0zE+M4m5EzOYOzGDsWkyqySWSViLKCdhfaZUVeVgXRuv7atj67569jtaAchPS+Sa6Tmcn2fjvFFWtHJrspgiYS2inIT12ap2edi2v57XD9RT6+7kSGMHFmMcl4xPZ/aEdGaOSyNZPqCMehLWIspJWA+mFq+fdyqcvHGggTcPNuDy+BmXYcJi1PPNcWnMGp/GRLtZVgaMQhLWIspJWA+V7p4PKHdVuthUXhtul6QnGfjmuDTmTszggjE2LEZ9hCsVIGEtop6E9XBpaO1kx2dO3jro5O0KJwVpJvZUu5mSY2FmYSrfGJdGSa6FOK2sWRIJEtYiyklYR0KgO8ieajdvf9bIOxVOyqvcBFUwGXRcXpxJcZaZiwtSKUg3SctkmEhYiygnYR0NWrx+3vu8kZ2HGnn7s0aONnuAUMvk4oJUZo1P4/w8G1mWhAhXOnJJWIsoJ2EdjY41edjZE97vfd6ELxCkrSvA6BQjX8tP4cKxoa8Ms1yUM1gkrEWUk7COdsFgkIP17bz3eRPvft7EB0eaaOsMkJdiRKtRmNET3BeOsZE+QsJ73bp1rFixggMHDlBWVtbvfRjz8vJISkpCq9Wi0+nCYdvc3MyiRYuorKwkLy+PF154AavVespzSliLKCdhHWu6gyp7a9x8cLiZ94808+GRZtq6AgBcMj6NDHM8F4yxccEYGzlWY4Sr/WoOHDiARqNh2bJlp7xpbl5eHrt27SI1NfWk7T//+c+x2WwsX76cVatW4XK5eOSRR055TglrEeVi84a55zKtRmFqrpWpuVaWfjOf7qDK/tpW3j/cxO4qF5s/cbD2wyoAspLjubQog3GZSVyQZ6Mw3YQmBq6sLCoqOqvHb9q0iR07dgCwZMkSZs2addqwFmKkk7COMK1GYXJOMpNzkoHQAlQH69v4sLKZDyub2XagnufePwqAOT4022R0SiKlo61MzbUQH6eNZPlnRVEULrvsMhRFYdmyZSxduhSA+vr68E1y7XY7DQ0NkSxTiKggYR1lNBqFIruZIruZGy7KQ1VVqpq9lFU2s6uymdZOP/+15SAAcVqFSVnJXDohnYJ0E9NGW4ftQ8s5c+ZQV1fXa/vKlStZuHDhgJ5j586dZGVl0dDQwNy5c5kwYQIzZ84ccA2rV69m9erVADidzgE/TohYJD3rGNTc4eNfR13sOurio6PNKCiUVTYDkG1J4IrJmWQlJzBttJWJdjN6XWQu1Jk1a9Ype9YnWrFiBSaTiXvvvZfx48ezY8cO7HY7DoeDWbNmcfDgwVM+XnrWIspJz/pcZEvUM2diBnMmZgDgCwTZV9vCv465+ddRF5/Vt/HUO0cA0Os0zB6fTo41gZJRFkpyLWRbEiJ+sU5HRwfBYJCkpCQ6OjrYunUrDz74IAALFixgzZo1LF++nDVr1gx4pC7ESCYj6xHK0eJl9zE3u4+5aO7w8/LHtXQFggCkmgxcWZxBujmeqbkWpuRYBnVlwQ0bNnDHHXfgdDqxWCyUlJSwZcsWamtrufnmm9m8eTOHDx/mqquuAiAQCPD973+f+++/H4CmpiauvfZajh07xqhRo1i3bh02m+2U55SRtYhyMnVPDIy/O8injjbKq1zsrnLT0RVgy7768P6xqYlcWpROliWBKTkWJmWZY+rDSwlrEeUkrMVX1+L1s7emhfIqN3uq3DhavHxSE1pdUKdR+Oa4NFJNBqbkJjMl28L4zKSI9b9PR8JaRDkJazG46lo62VPtZk+Vi4N17Xx0zIXb4wdg+igrvu4gxdnJTM4OTTccl2HCoIv8CFzCWkQ5CWsxtFRVpdrl5ZOaFiqbOth5qJFPqlto7QxwQZ6N3VUuxmUkUZyVzJScZCbYzUy0m0nQD2+AS1iLKCdhLYbf8bnfn9a1srvKzd6aFqpcHiobQ6sNahSYNT4Nc3wck7KSmZhlZqI9CWuiYchqkrAWUU7CWkQHVVVxtHSyt6aFvbWtuDw+tu2vx9HSCYRuSNzpD1JkNzMxy8yUnGTGZySRYx2caYQS1iLKSViL6NbU3sV+RysHHK3sqw19HXa2Y7ckUOPykmTQ9VyxmcTELDPjM82Mz0g64zaKhLWIchLWIvZ4fd18WtfKp3Vt7K8NBXlti5dad2gUrigwJiWRIruZCZlJFGcnU5BuOuUoXMJaRDm5glHEngS9lvNGWTlv1BdrVHd3B6l2ezngaOOAo5VP61r5pKaFVz5xUJhhoqK+HZNBx7gMExN6Qnx8RhLjM5PkpsXinCAjaxHV2rsCHKxr5WBdOwfrWjlQ18bBujZavKHphDPG2Dja5OHQ6p9w9/+8yLiMJCZkJpGflkiCXsYiImrIyFqMbCaDjumjbUwf/cXl5qqqUt/axYG6Vo41edhT5eZAUOXZdyvx9VxSP2OMjYa2LgrTTYzPTKIwIzQSH5OaGLUX9ghxKjKyFiNCaWkp739QRmWTh8/q26hq9lBe5eaz+jYqmzx0B0Nv3dLRVtxeP+MyTBSmJ1GYYWJcehJ5EuJiaMnIWojjdFoNBekmCtJNJ23vCnRz2NnBZ/Vt1LV0suuoi/21rby6tw5VhUS9ls5AkLwUI+MykihMN1FkNzM6JZGxaYkxtUaKGLlkZC1GhK8yG6TT382hhnYONbRT0dBGRX3o+8qmDs7Ps/HBkWYUBUbZjBSkmchPT6QwPYn8nv8hmOMHb6VCMeLJyFqIryo+TktxdjLF2cknbe/0d1PZ2EFFT5AfcrZzqL6dg/VtVLuOhI9LSzJQkBYK7gn2JEbbEslPTyTTHB/x9cLFyCMjazEiDMc860B3kCqXl8+PB3hDO5872znsbMfj68bfHfrrYdRrmdtzY4j8NBNj0xLJTzORl2KUGSrnLrkoRgiI7EUxwWCQxnYfh5ztfO7s4POGdjq6Arz7eRM1bm/4uBljbFS7vOHwPv7nmNTQaDwW7lwvvjJpgwgRaRqNhnRzPOnmeL6Wn3rSPo8vwJHGDg47O2ho62RPVQuHG9t5YVcV5ngdda1dACTEaclLTeT8PCvJCXGMSU1kTGoiY1NNJBulNy4krMUgW7duHStWrODAgQOUlZX1ebPcgwcPsmjRovDPhw8f5le/+hU/+9nPWLFiBX/6059IS0sD4OGHH+bKK68ctvoHm1GvY1JWMpOyTu6Lq6pKXWsnR5wdHG7s6An0dhpau/jrB8fCUw0hdM/NC/JsJBp0jE1LJC8lFOSjUxJINEiQnyskrMWgKi4u5sUXX2TZsmX9HjN+/HjKy8sB6O7uJjs7O3w/RoC77rqLe++9d8hrjSRFUbAnJ2BPTuBrBSePxn2BIFUuT0+Qt3OksYMadye7q1z8v39Vh48rHW2lyuUJh3deaijIx6YlMspmlCmHI4yEtRhURUVFZ3T89u3byc/PZ/To0UNUUezR6zTkp5nITzMBGSft6+gKUNnUQWWjB2d7J59Ut1LZ1MHr++tp6vABod74B0easSfHMzrF2DMKDwX6KJuR0SlGjPJBZ8yR/2IiotauXct111130rbHH3+cv/zlL5SWlvLb3/4Wq9Xa52NXr17N6tWrAXA6nUNeazRINPTdVgFo7fRT2dhBjcvLxQWpVDZ1cLTJw9Z99bR4fQRVON5dSU8yMCnLTIrJQF6KkdEpiYxOMTLKZpSFsaKUzAYRZ2zOnDnU1dX12r5y5UoWLlwIwKxZs3j00Uf77Fkf5/P5yMrKYt++fWRkhEaQ9fX1pKamoigKDzzwAA6Hg2eeeea0NckSqafW6vVzrNkTDvDKxg68vm7KKptpaOsKHzdjjI2D9W2MthkZlZLIaJuR/LREsiwJjE5JJD3JILNWvhqZDSKG37Zt2wbleV599VWmTZsWDmrgpO9vueUW5s+fPyjnOteZE+L6vAAIQjNWjjWHbsvW2NZJfp2JY00eyqtcvPJxLVNzktld1QKAQach12ZktM1IQbqJzOR4RtlCI/Jc6ZMPKQlrETHPP/98rxaIw+HAbrcDsGHDBoqLiyNR2jnFqNcxIdPMhExzr33+7iA1Li9Hmz0ca/ZwrKmDY80ejjZ5cLZ38XF1y0nHZ5gNFKYnkW42hALcamRUT3slzaRHo5HFsr4qaYOIQbVhwwbuuOMOnE4nFouFkpIStmzZQm1tLTfffDObN28GwOPxkJuby+HDh0lO/mK098Mf/pDy8nIURSEvL48//vGP4fA+FWmDDD9VVWnq8HGs2UNVT4Afa/bQ3hVgT5WbutZOToyXRIOWTHM8ucdH4tZQiyUjObRthK+1IlcwCgES1tGoK9BNjcvLsWbPF6PzJg9VrlCot3UGmJprYU+VG4DkhDhybQmcl2slQa8lx5pArtVIjjWBHKvxjO/LGWWkZy2EiE4GnZaxaSbGppn63N/iCX3oWeXyUO3yUNXspcrlwdnWxZsHG+jquZHEcTPG2OgKBMPhHfrzi+9Her9cwloIERHJxjgmG5OZnNP7Q09VVXG2d1HV7KXa5aHa5cXTFWBPdQt7a1rYsq8uvHAWgFaB1CQD9uSTAzwvxUhmcjzZlpgfmUtYCyGij6IopCfFk54Uz/TRvefZB4MqDW1doRG5y0Odu5PKJg81bu9JYX5Bno2yymYgdNl+jjWBbEs8OdZEsizxZFtCwZ5tTSA5Ibp75hLWQoiYo9EoZCbHk5kcT2merdf+42Fe6/b2tFm8VLu81Li91Li8vPGps1eb5YIxNlq9frIsCWRbEsi2hv7MsiSQbY0n3RTZlRElrIUQI86JYT6tj5H58ZksNScEuMcX4JOaVmrcXj466qLF6w8fP320lY+r3diTE8iyxIcDPfxncjxZ1oQhvYxfwloIcc5RFIVUk4FUk4GpuZY+j2nr9FPr7qTW7cXZ3sUFY2zU9gT7+583UdfaSVAFk15Lu68bAIsxjqzk4yEeCvULx6b0e44zIWEtRoTU1NTTHyTEGUiKj2N8ZhzjM5P63B/oDlLX2onD7aXG3UmN24ujxUutu5Nql4cPjjTR1hng7rnjBiWsZZ61EEIMkbZOPypgjo+Ti2KEECIGnHVYy4X6QggRAySshRAiBkhYCyFEDJCwFkKIGCBhLYQQMUDCWgghYoCEtRBCxIAzvYJR7pQphBARICNrIYSIARLWQggRAySshRAiBkhYCyFEDJCwFkKIGCBhLYQQMUDCWgghYoCEtRBCxAAJayGEiAES1kIIEQP+fwohKgDxYSbCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<sympy.plotting.plot.Plot at 0x11b5e0978>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.plot(*[sympy.re(r) for r in roots], (K_c, 1, 20))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that one root gets a positive real part around $K_c=12.5$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using the control library\n",
    "We quickly run out of SymPy's abilities when calculating closed loop transfer functions. Let's try to do it with the control library instead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import control\n",
    "import numpy\n",
    "import scipy.signal\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = control.tf([1, 0], [1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "G_v = 1/(2*s + 1) \n",
    "G_p = G_d = 1/(5*s + 1)\n",
    "G_m = 1/(s + 1)\n",
    "K_m = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tbcontrol.loops import feedback"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def closedloop(K_c):\n",
    "    G_c = K_c\n",
    "\n",
    "    G_CL = K_m*feedback(G_c*G_v*G_p, G_m)\n",
    "    return G_CL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\frac{20 s^3 + 34 s^2 + 16 s + 2}{100 s^5 + 240 s^4 + 209 s^3 + 103 s^2 + 29 s + 3}$$"
      ],
      "text/plain": [
       "\n",
       "           20 s^3 + 34 s^2 + 16 s + 2\n",
       "------------------------------------------------\n",
       "100 s^5 + 240 s^4 + 209 s^3 + 103 s^2 + 29 s + 3"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "closedloop(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipywidgets import interact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "smootht = numpy.linspace(0, 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "loop = G_v*G_p*G_m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alchemyst/anaconda3/lib/python3.7/site-packages/matplotlib/figure.py:98: MatplotlibDeprecationWarning: \n",
      "Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.\n",
      "  \"Adding an axes using the same arguments as a previous axes \"\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmUlOWd9vHvrzd2mqUbEUFARUUUUXqMYo7RiIo6GqMhAQ1gjKKJRuOSTHIyb5xkJm8yrwuaaFREx7giigsxGoia6ERcaBQQIyoQxBaUbpaGZu3l9/5RD0WDTVMFXXXXcn3OqdPPXfV01UWFcPls92PujoiISKIKQgcQEZHsouIQEZGkqDhERCQpKg4REUmKikNERJKi4hARkaSoOEREJCkqDskrZrbMzDabWZ2ZfWZmD5hZ5zZ43wFm5mZW1Mo6/2FmD+/rZ4mEpuKQfHSOu3cGhgHHAD8NnEckq6g4JG+5+2fATGIFAoCZlZrZg2ZWbWYfm9m/m1lB9FpBNP7YzFZF65VGv/pq9HNdtDVzQjJZzGywmf3NzNaZ2Xtmdm6z1zqY2S3R59aa2d+j5042s6pd3meZmY2Mlo8zs0ozW29mn5vZrcl/SyJfpOKQvGVmfYEzgcXNnv4dUAocBHwFGA98J3rt4uhxSvR6Z+CO6LWTop/d3L2zu7+eRI5i4I/ALKAX8APgETM7LFrlZmA4MALoAfwYaErgrW8Hbnf3rsDBwLREM4m0RsUh+egZM9sAfAKsAm4EMLNC4FvAT919g7svA24BxkW/dxFwq7svdfc6Yru4xrR2XCNBxxMrod+4+zZ3fxl4Dhgbbe1cAlzj7p+6e6O7z3b3rQm8bz1wiJmVuXudu7+xjzlFABWH5Kfz3L0LcDJwOFAWPV8GlAAfN1v3Y+CAaLlPC68VAfvtY54+wCfu3nwrYvvnlgHtgSV78b7fBQ4FFpnZHDP7133MKQKoOCSPufsrwAPEdgUB1BD7r/T+zVY7EPg0Wl7RwmsNwOfAvkwzvQLot/1Yyi6fWwNsIbaraVcbgY7bB9EWU/n2sbt/5O5jie3++m/gSTPrtA85RQAVh8htwGlmNszdG4kdB/iVmXUxs/7AdcD2U2gfA641s4HRKbz/F3jc3RuAamLHHQ7aw+cVmFn7Zo92wJvESuDHZlZsZicD5wBTo62Q+4FbzayPmRWa2QnR730ItDezs6PjJP8OtNv+QWb2bTMrj95jXfR04z59WyKoOCTPuXs18CDwf6KnfkDsH/GlwN+BR4n9w0308yFiZ1D9k9iWwA+i99kE/Ap4LToz6vjdfORYYHOzxxJ33wacS+xAfQ3we2C8uy+KfucG4F1gDrCG2NZDgbvXAt8HphDbOtkIND/LahTwnpnVETtQPsbdtyT5FYl8gelGTiIikgxtcYiISFJUHCIikhQVh4iIJEXFISIiSdnXK14zUllZmQ8YMCB0DBFJwty5cxk+fHhK3ru+qorGdbUAWIf2tDu4pcti8tvcuXNr3L18z2sGPqvKzO4H/hVY5e5HtvD6ycCzxE59BHjK3X+5p/etqKjwysrKtowqIilmZqTi36PN8+axbMzY+Lj/ww/RsaKizT8n25nZXHdP6IsJvcXxALFJ4h5sZZ3/dXdNlSAiSXN3Pvv1r+PjLmecodJoA0GPcbj7q8QuaBIRaXPrn/sTW+YvAMCKi+l1w/WBE+WGbDg4foKZzTezF8xsyO5WMrOJ0b0HKqurq9OZT0QyUNOWLayatOMWJD0mjKekX7+AiXJHphfH20B/dz+a2H0Sntndiu4+2d0r3L2ivDyh4zsiksPWPPQQDStWAlDYowc9L788cKLckdHF4e7ro/se4O7PA8VmVraHXxORPNewZg2r75kcH5dddSWFXboETJRbMro4zKy3mVm0fByxvKvDphKRTFdzx5001dUBUDJwIN1Hjw6cKLcEPavKzB4jdjOdsujeyTcCxQDufjfwDeB7ZtZAbCbRMa5ZGUWkFVuX/pO1jz8eH/f60Y+w4uKAiXJP0OKIbjLT2ut3sOOeziIie7Tq5puhMXbbkY7HHUfnU04OGygHZfSuKhGRZGx86y3qXn45Pu714x8T7e2WNqTiEJGc4E1NrPp/N8XHXc89hw5H7vYMftkHKg4RyQnr//Q8WxYuBMDataPXD38YOFHuUnGISNZr2rqV6kmT4uMe48dT3KdPwES5TcUhIllv7cMPU79iBQCF3bvTc+JlgRPlNhWHiGS1hrVrqbn7nvi47Epd7JdqKg4RyWo1d91F04YNAJQMGED3b30zcKLcp+IQkay1bfly1j42NT4uv/46XeyXBioOEclaq26dBPX1AHQYPpwuI0cGTpQfVBwikpU2z5vHhj//OT7e78c/0sV+aaLiEJGs4+58ftPN8XGXUaPocPTRARPlFxWHiGSdupdeYvPcubFBcTG9rrs2bKA8o+IQkazi9fWsuvmW+Lj72DGUHHhgwET5R8UhIlll3fTpbFu2DICCzp0p+973wgbKQyoOEckajXUbqf7djjst9Lx8IkXduwdMlJ9UHCKSNdbcfz+Nq2M3AS3q3Zse48YFTpSfVBwikhUaqqtZ/cAD8XH5NddQ0L59uEB5TMUhIlmh+o478U2bAGh32GGUnntO4ET5S8UhIhlv69KlrHvyyfi41w3XY4WFARPlNxWHiGS8VbfeuuM+4iccT6cvfzlwovym4hCRjLbp7Xeoe/Gl+LjX9TdoapHAVBwikrHcnVU375hapOvZZ+s+4hlAxSEiGavu5ZfZ/PbbsUFxMeXX6j7imUDFISIZyRsaYtOmR7qPGUNJ374BE8l2Kg4RyUi1zzzDtiVLACjo1Imy710ROJFsp+IQkYzTtHnzzlOLXHYpRT16BEwkzQUtDjO738xWmdnC3bxuZvZbM1tsZgvM7Nh0ZxSR9Fvz8MM0fP45AIXlZfQYPz5wImku9BbHA8CoVl4/ExgUPSYCd6Uhk4gE1Fhby+p7p8TH5VdeSUHHjgETya6CFoe7vwqsaWWVrwEPeswbQDcz2z896UQkhNX33kvT+vUAlPTvT7cLLgicSHYVeotjTw4APmk2roqe+wIzm2hmlWZWWV1dnZZwItK26j/7jDUPPRwfl//wGqy4OGAiaUmmF0dLl4d6Syu6+2R3r3D3ivLy8hTHEpFUqLnz9/jWrQC0HzKELmecETiRtCTTi6MK6Nds3BdYESiLiKTQgOIS1j31VHzc6/rrsIJM/ycqP2X6/yozgPHR2VXHA7XuvjJ0KBFpe9eUl+08keGIEYETye4UhfxwM3sMOBkoM7Mq4EagGMDd7waeB84CFgObgO+ESSoiqbT53YWc0aVrfNzruusCppE9CVoc7j52D687cGWa4ohIINWTdkwt0uX00+lw1FEB08ieZPquKhHJcRvfeJONs2fHBgUFlP/wmrCBZI9UHCISjLvvtLVR+vXzaHfQQQETSSJUHCISTN1f/8bm+fMB2NbURPmV2jOdDVQcIhKENzVRffvt8fHjteso7tMnYCJJlIpDRIJY/8ILbP3gAwCsY0cmr14dOJEkSsUhImnnDQ3UNJs2vce4cayOruGQzKfiEJG0q312BtuWLQOgoEsXel6iS7SyiYpDRNLKt22j5s474+Oel3yHwtLSgIkkWSoOEUmrddOnU78iNuVcYffudB+nmzRlGxWHiKRN09at1Nx9T3zc89JLKezcKWAi2RsqDhFJm3WPT9vplrDdL2x11iHJUCoOEUmLpi1bqLl3cnxcdtlECjp0CJhI9paKQ0TSYu3UqTRW1wBQtN9+dPvWNwMnkr2l4hCRlGvavJnVU+6Lj3tOvIyCdu0CJpJ9oeIQkZRbO/VxGmuirY3evek2enTgRLIvVBwiklKxrY0p8XHZ5RMpKCkJmEj2lYpDRFJq3bRpNEbzUBX17k3pBRcETiT7SsUhIinTtGULNc22Nnpedqm2NnKAikNEUmbd9Ok7zqTq1Ytu3/hG4ETSFlQcIpISvm3bzmdSXXqpzqTKESoOEUmJdc8+S8PKlQAU9uxJt9Ha2sgVKg4RaXPe0MDqe5sd2/jOxbpKPIeoOESkza2fOZP65csBKCgtpdsYzUmVS1QcItKm3H2nrY0eF12kGXBzjIpDRNrUxv/9X7YuWgSAdehA93HfDpxI2pqKQ0Ta1OrJ98aXu43+BkXduwdMI6kQtDjMbJSZfWBmi83sJy28frGZVZvZvOhxaYicIpKYzQsWsKmyMjYoKqLnxRcHzSOpURTqg82sELgTOA2oAuaY2Qx3/8cuqz7u7lelPaCIJG3NA3+IL5eefRbFffoETCOpEnKL4zhgsbsvdfdtwFTgawHziMg+qF+xgvUzZ8bHPbS1kbNCFscBwCfNxlXRc7u6wMwWmNmTZtZvd29mZhPNrNLMKqurq9s6q4jswZqHHobGRgA6Hn887QcPDpxIUiVkcVgLz/ku4z8CA9x9KPAi8Icv/kr0i+6T3b3C3SvKy8vbMKaI7EljXR3rnngiPu5x8YSAaSTVQhZHFdB8C6IvsKL5Cu6+2t23RsN7geFpyiYiSaidPp2mujoASgYOpPNJJwVOJKkUsjjmAIPMbKCZlQBjgBnNVzCz/ZsNzwXeT2M+EUmAu7P2sanxcY8JE7ACnemfy4KdVeXuDWZ2FTATKATud/f3zOyXQKW7zwCuNrNzgQZgDXBxqLwi0rJNb81h27JlABR07kzpueeEDSQpF6w4ANz9eeD5XZ77ebPlnwI/TXcuEUncumnT4sul555DQceOAdNIOmh7UkT2WsPatWyYNSs+7vbNbwZMI+mi4hCRvVb79DN4fT0A7Y8eSvvDDw+cSNJBxSEie8Xdd9pN1V1bG3lDxSEie2XXg+JdzzwzbCBJGxWHiOyVdY8/Hl/WQfH8ouIQkaQ1rlvHhr/8JT7WQfH8klBxRDPZiogAsH7mrB0HxY88UgfF80yiWxyLzewmMzsipWlEJCusf+65+LIu+Ms/iRbHUOBDYIqZvRHNRNs1hblEJEPVr1q142ZNBQU6KJ6HEioOd9/g7ve6+wjgx8CNwEoz+4OZHZLShCKSUepefhk8NpF1x3/5F4o0G3XeSfgYh5mda2ZPA7cDtwAHEZv2/PlWf1lEcsqGWTsOincZOTJgEgkl0bmqPgL+Ctzk7rObPf+kmWn+ZJE80bh+PRvfeis+7jLy1IBpJJQ9Fkd0RtUD7v7Lll5396vbPJWIZKSNb74JDQ0AtB8yhOL999/Db0gu2uOuKndvBE5JQxYRyXCbXn8jvtzpxBMDJpGQEt1VNdvM7gAeBzZuf9Ld305JKhHJSBtffz2+3GnECQGTSEiJFseI6Gfz3VUOfLVt44hIpqpfuZJt//wnANauHR2OOSZwIgkloeJwd+2qEslzm5odFO84fDgF7doFTCMhJXwHQDM7GxgCtN/+3O4OmItI7tk8f358uUPF8IBJJLREr+O4G/gW8APAgNFA/xTmEpEMs/ndhfHlDkOPDphEQkt0ypER7j4eWOvuvwBOAPqlLpaIZBJvaGDrhx/Gx+2HaNq6fJZocWyOfm4ysz5APTAwNZFEJNPUV1XhW7cCUNSrF0XduwdOJCEleozjOTPrBtwEvE3sjKopKUslIhlla3Q2FUDJQP03Y75L9Kyq/4wWp5vZc0B7d69NXSwRyST1y5fHl0sGDAgXRDJCMmdVjQAGbP8dM8PdH0xRLhHJIA3r1sWXNRuuJFQcZvYQcDAwD2iMnnZAxSGSB5pq18eXC7t2CZhEMkGiWxwVwBHu0ST8IpJXGjdsiC8XdNU93PJdomdVLQR6t/WHm9koM/vAzBab2U9aeL2dmT0evf6mmQ1o6wwi6XD3K0uYvaRmp+dmL6nh7leWZEWWpvXNtzi+WByZ9OeT1Et0i6MM+IeZvQVs3f6ku5+7tx8cTdd+J3AaUAXMMbMZ7v6PZqt9l9i1I4eY2Rjgv4ldiCiSVYb2LeWqR9/hjguPYcTBZcxeUhMfh8jy/elPMvqEQgaWdeKfNRuZVvkJ36zox6Pvv9ni7wz6fCmdo+WNs1+nfsXKnV7/l5o6pt1bRbuKvgws68xH3fty1dvbgvz5JPUskb1PZvaVlp5391f2+oPNTgD+w93PiMY/jd7z183WmRmt87qZFQGfAeV72mVW0afQKyd2bm0VkeBmj1vCiIPLgnz21bN+zl9XPp3w+rfc20C/mj2vt91TR43ilN/8LKk/X3TCTeIfIm3KzOa6e0Ui6yZ6Ou5eF0QrDgA+aTauAr60u3XcvcHMaoGewBf+CpvZRGAiwPD9E90DJxLOiYeEOztp/4v2p+dpPRNef21no19N4v+or10yf6/+fGaW9O9I+rVaHGb2d3f/spltIHYWVfwlwN19X46StfQ3ZNe/mYmsE3vSfTIwGWJbHPuQSyQtXltcHWyL43dvPsmUOS9y6H5d+PDzDZx0aDm9u7bf7frTzp/Fkter6bAVRvYfSXmHL5bCyvVbePXDag7brwtrCvvx2hPTtMWRRZIp7VaLw92/HP1Mxfl3Vew831VfYMVu1qmKdlWVAmv2+M59joH/qGyjmCL7rvkxjV2PcaS7PGYvqeF/ZpVy14X/lXCWK+o+4bF2sf/rDT/1GxzV96QvvOdVj77DHZNi79E14J9PUi/R2XF7tPAo3sfPngMMMrOBZlYCjAFm7LLODGBCtPwN4GWdEizZaEFV7U7/iI44uIw7LjyGBVXpn4Bhb7J0Ldmxc6F26xfXy6Q/n6ReomdVvU3sv/zXEtt91A1YaWargMvcfW6yHxwds7gKmAkUAve7+3tm9kug0t1nAPcBD5nZYmJbGmOS/RyRTHDFVw7+wnMjDi4L8l/je5OleXFs2LbhC69n0p9PUi/R4vgz8LS7zwQws9OBUcA04Pd88aB2Qtz9eeD5XZ77ebPlLcTu/SEiAXUp2bG3ev229a2sKfkg0dOPKraXBoC7zwJOcvc3AN0/UiTHlZaUxpfXbNnzYUbJbYlucawxs38DpkbjbwFro4v4mlKSTEQyRt8ufePLyzcsb2VNyQeJbnFcSOysp2eAZ4EDo+cKgW+mJpqIZIoBXQfEl5fVLguWQzJDohcA1hC733hLFrddHBHJRAd2PZAiK6LBG/i07lPqttXRuUSzM+SrRE/HLTezm8zseTN7efsj1eFEJDOUFJYwsNuOO/8tWrMoYBoJLdFdVY8Ai4jdZ/wXwDJi12GISJ4Y0nNIfPm91e8FTCKhJVocPd39PqDe3V9x90uA41OYS0QyzNDyofHluZ8nfemW5JBEz6qqj36uNLOziU0N0reV9UUkx1Tst2Pi1DmfzaGhqYGigoTvPi05JNEtjv8ys1LgeuAGYApwbcpSiUjGGdB1APt13A+Auvo67a7KYwkVh7s/5+617r7Q3U9x9+HRlCAikifMjOP337GH+o0VbwRMIyElelbVQDO71cyeMrMZ2x+pDicimeX4PjuK4/WVrwdMIiEluoPyGWITDv4RXSkukreab3HMWzWPdVvW0a19t4CJJIREj3Fscfffuvtfo7OqXknRXQFFJIOVdSiLn13V6I38repvYQNJEIkWx+1mdqOZnWBmx25/pDSZiGSkUw88Nb780vKXAiaRUBLdVXUUMA74Kjt2VXk0FpE8MvLAkUyaOwmA2Z/OZsO2DTtNuy65L9Etjq8DB7n7V6Kzqk5xd5WGSB46sOuBDO4xGIBtTdu01ZGHEi2O+cTu+iciwlkDz4ovP7f0uYBJJIREi2M/YJGZzdTpuCIyauAoDAPgrZVvUbWhKnAiSadEj3HcmNIUIpJVenfqzYgDRvDap6/hOE999BRXH3t16FiSJoleOf5KS49UhxORzDV60Oj48tOLn6a+qb6VtSWXtFocZrbBzNa38NhgZrpjvUgeO6nfSZR3KAegZnMNr3yi/5bMF60Wh7t3cfeuLTy6uHvXdIUUkcxTXFDMeYecFx8/+eGTAdNIOiV6cFxE5AsuOPSC+EHy2Stm6yB5nlBxiMheO6DzAYw4YARA/CC55D4Vh4jsk9GH6iB5vlFxiMg+OanvzgfJX17+cuBEkmpBisPMepjZX8zso+hn992s12hm86KHLjgUyUDFBcWcP+j8+PjBfzwYMI2kQ6gtjp8AL7n7IOClaNySze4+LHqcm754IpKMMYePobigGIAF1QuYt2pe4ESSSqGK42vAH6LlPwDntbKuiGS4sg5lnH3Q2fGxtjpyW6ji2M/dVwJEP3vtZr32ZlZpZm+YWavlYmYTo3Urq6ur2zqviOzBuCPGxZdfWv6STs3NYSkrDjN70cwWtvD4WhJvc6C7VwAXAreZ2cG7W9HdJ7t7hbtXlJeX73N+EUnOod0PZUSf2Km5Td7EI+8/EjiRpErKisPdR7r7kS08ngU+N7P9AaKfq3bzHiuin0uBvwHHpCqviOy78UeMjy9P/2g6tVtrA6aRVAm1q2oGMCFangA8u+sKZtbdzNpFy2XAicA/0pZQRJI2os8IDul2CACbGzbz2KLHAieSVAhVHL8BTjOzj4DTojFmVmFmU6J1BgOVZjYf+CvwG3dXcYhkMDPju0d9Nz5+5P1H2FS/KWAiSYUgxeHuq939VHcfFP1cEz1f6e6XRsuz3f0odz86+nlfiKwikpxRA0ZxQOcDAFi3dR1PL346cCJpa7pyXETaVFFBEROGTIiPH3jvAeobNQ1JLlFxiEibO++Q8+jRvgcAn238TPclzzEqDhFpcx2KOux0XceUd6fQ0NQQMJG0JRWHiKTEmMPG0LUkdr+35RuWM3PZzMCJpK2oOEQkJTqXdObbg78dH09eMJkmbwqYSNqKikNEUubCwRfSqbgTAEtrlzLr41mBE0lbUHGISMqUtivlwsMvjI/vmX+PtjpygIpDRFJq/BHj6VjUEYDF6xbzl4//EjiR7CsVh4ikVLf23Rh7+Nj4+O75d2urI8upOEQk5SYMmbDTVoeOdWQ3FYeIpFz39t25cPCOYx13z7ubxqbGgIlkX6g4RCQtJhwxIX6G1ZLaJbquI4upOEQkLbq178ZFgy+Kj++af5euJs9SKg4RSZvxR4ynS3EXAJatX8aflv4pcCLZGyoOEUmb0naljBuyYw6ru+bfRX2TZs7NNioOEUmrcYPHUdquFIBP6z7l6Y90v45so+IQkbTqXNKZS468JD6evGAyWxu3BkwkyVJxiEjajTlsDD3b9wTg802fM+2DaYETSTJUHCKSdh2LO3LZ0Mvi4ynvTqGgnf45yhb6X0pEghh96Gh6d+oNwJota+h5Ws/AiSRRKg4RCaKksITvH/39+LjsrDJqt9YGTCSJUnGISDDnHHwOA7oOAKCwYyH/s/B/wgaShKg4RCSYooIirjzmyvj4kfcfoXpTdcBEkggVh4gEdXr/0xncYzAAWxq3cM+CewInkj1RcYhIUAVWwNXHXh0fT/9wOsvXLw+YSPZExSEiwZ3Y50Q2LtoIQIM3cMc7dwROJK0JUhxmNtrM3jOzJjOraGW9UWb2gZktNrOfpDOjiKSPmfHZE5/Fxy8se4FFaxYFTCStCbXFsRA4H3h1dyuYWSFwJ3AmcAQw1syOSE88EUm3zUs2c0q/U+Lj296+LWAaaU2Q4nD39939gz2sdhyw2N2Xuvs2YCrwtdSnE5FQrj7magos9s/Sa5++xpzP5gROJC3J5GMcBwCfNBtXRc+1yMwmmlmlmVVWV+t0PpFsdEj3QzjnoHPi49vm3oa7B0wkLUlZcZjZi2a2sIVHolsN1sJzu/0b5O6T3b3C3SvKy8v3LrSIBHflsCspKSgBYEHNAl5a/lLgRLKrlBWHu4909yNbeDyb4FtUAf2ajfsCK9o+qYhkkv0778+Yw8fEx7e/fbtuMZthMnlX1RxgkJkNNLMSYAwwI3AmEUmDy466jM7FnYHYLWafXqybPWWSUKfjft3MqoATgD+Z2czo+T5m9jyAuzcAVwEzgfeBae7+Xoi8IpJe3dp347tHfTc+vmveXWxu2BwwkTQX6qyqp929r7u3c/f93P2M6PkV7n5Ws/Wed/dD3f1gd/9ViKwiEsZFgy+ivEPseGX15moe/sfDgRPJdpm8q0pE8liHog58f9iOadfvX3g/a7esDZhItlNxiEjGOu+Q8xhYOhCAuvo6Ji+YHDiRgIpDRDJYUUER1xx7TXw89YOpfLLhk1Z+Q9JBxSEiGe2r/b7KMb2OAaChqYHfvfO7wIlExSEiGc3MuG74dfHxC/98gfdW6wTLkFQcIpLxhvUaxqkHnhofT6qcpKlIAlJxiEhWuObYayi0QgDe/OxN/v7p3wMnyl8qDhHJCgNLB3LBoAvi40lvT6KxqTFgovyl4hCRrPG9Yd+jQ1EHAD5a+xEzlmgWohBUHCKSNco6lPGdId+Jj++Yd4emIglAxSEiWWXCkAn0bN8TgFWbVmkqkgBUHCKSVToWd9xpKpL7Ft7Hmi1rAibKPyoOEck65w86n4NKDwJgY/1G7p5/d+BE+UXFISJZp6igiGuHXxsfP/HBE3y8/uOAifKLikNEstJX+n6Fiv0qAGjwBm6be1vgRPlDxSEiWcnMuKHihvj4xeUv8s6qdwImyh8qDhHJWkPKhnDmgDPj45srb9ZUJGmg4hCRrHb1sVdTXFAMwILqBcz6eFbgRLlPxSEiWa1vl75cePiF8fFtc29jW+O2gIlyn4pDRLLeZUMvo2tJVwCq6qqYumhq4ES5TcUhIlmvtF0pVxx9RXw8+d3J1G6tDZgot6k4RCQnjDlsDH079wWgdmst9y64N3Ci3KXiEJGcUFxYzA+H/zA+fnTRo7o/eYqoOEQkZ5ze/3SGlg8FoL6pnt++/dvAiXKTikNEcoaZ8aOKH8XHf172Z+ZXzw+YKDepOEQkpwzrNYzT+p8WH99SeYsuCmxjQYrDzEab2Xtm1mRmFa2st8zM3jWzeWZWmc6MIpK9rj32WooKigB4Z9U7vLj8xcCJckuoLY6FwPnAqwmse4q7D3P33RaMiEhz/br2Y+zhY+PjSXMnUd9YHzBRbglSHO7+vrt/EOKzRSQ/XD708vhFgZ9s+IRHFz0aOFHuyPRjHA7MMrO5ZjaxtRXNbKKZVZpZZXV1dZriiUimKm1XyuVDL4+dqDvPAAAEdUlEQVSP75l/D2u3rA2YKHekrDjM7EUzW9jC42tJvM2J7n4scCZwpZmdtLsV3X2yu1e4e0V5efk+5xeR7Df28LH079ofgA31G/j9vN8HTpQbUlYc7j7S3Y9s4fFsEu+xIvq5CngaOC5VeUUk9xQXFnP98Ovj4yc+fIIl65YETJQbMnZXlZl1MrMu25eB04kdVBcRSdjJ/U7mS/t/CYBGb+SmypsCJ8p+FuL8ZjP7OvA7oBxYB8xz9zPMrA8wxd3PMrODiG1lABQBj7r7rxJ8/2ogXTcgLgNq0vRZ2UbfTcv0vbRM30vL0vW99Hf3hPbzBymOXGJmlTpVuGX6blqm76Vl+l5alonfS8buqhIRkcyk4hARkaSoOPbd5NABMpi+m5bpe2mZvpeWZdz3omMcIiKSFG1xiIhIUlQcIiKSFBVHGzCz/zSzBdH077Oi61HynpndZGaLou/maTPrFjpTpkj01gL5wsxGmdkHZrbYzH4SOk8mMLP7zWyVmWXchc8qjrZxk7sPdfdhwHPAz0MHyhB/AY5096HAh8BPA+fJJMncWiCnmVkhcCexOemOAMaa2RFhU2WEB4BRoUO0RMXRBtx9fbNhJ2Kz+uY9d5/l7g3R8A2gb8g8mUS3FtjJccBid1/q7tuAqUAyk6HmJHd/FVgTOkdLikIHyBVm9itgPFALnBI4Tia6BHg8dAjJSAcAnzQbVwFfCpRFEqDiSJCZvQj0buGln7n7s+7+M+BnZvZT4CrgxrQGDGRP30u0zs+ABuCRdGYLLZHvRgCwFp7TVnsGU3EkyN1HJrjqo8CfyJPi2NP3YmYTgH8FTvU8u2goib8z+a4K6Nds3BdYESiLJEDHONqAmQ1qNjwXWBQqSyYxs1HAvwHnuvum0HkkY80BBpnZQDMrAcYAMwJnklboyvE2YGbTgcOAJmLTuV/h7p+GTRWemS0G2gGro6fecPcrAkbKGLu7tUDYVOGY2VnAbUAhcH+it1DIZWb2GHAysWnVPwdudPf7goaKqDhERCQp2lUlIiJJUXGIiEhSVBwiIpIUFYeIiCRFxSEiIklRcYi0ATNrjGZHXmhmf9yXmYDNbJmZlbVlPpG2pOIQaRub3X2Yux9JbGK6K0MHEkkVFYdI23ud2MR9AJjZj8xsTnRfkl80e/4ZM5sb3ZdjYpCkIntBxSHShqJ7S5xKNGWGmZ0ODCI2dfgwYLiZnRStfom7DwcqgKvNrGeAyCJJU3GItI0OZjaP2PQqPYjdxArg9OjxDvA2cDixIoFYWcwndq+Sfs2eF8loKg6RtrE5ugNkf6CEHcc4DPh1dPxjmLsf4u73mdnJwEjgBHc/mlixtA8RXCRZKg6RNuTutcDVwA1mVgzMBC4xs84AZnaAmfUCSoG17r7JzA4Hjg8WWiRJuh+HSBtz93eiXVBj3P0hMxsMvG5mAHXAt4E/A1eY2QLgA2K7q0SygmbHFRGRpGhXlYiIJEXFISIiSVFxiIhIUlQcIiKSFBWHiIgkRcUhIiJJUXGIiEhS/j8a2Y6P+uLKegAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = control.rlocus(loop)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def response(K_C):\n",
    "    G_CL = closedloop(K_C)\n",
    "    poles = G_CL.pole()\n",
    "    plt.plot(*control.step_response(G_CL, smootht))\n",
    "    _ = control.rlocus(loop)\n",
    "    plt.scatter(poles.real, poles.imag, color='red')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5f3e593d9f7a4ff58ef56059318652f7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=10.5, description='K_C', max=20.0, min=1.0), Output()), _dom_classes=(…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.response(K_C)>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "interact(response, K_C=(1., 20.))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the step response is calculated quickly enough that we can interact with the graphics using the slider!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Direct substitution\n",
    "\n",
    "From our exploration above it is clear there is a zero of the characteristic equation at $K_C \\approx 13$. Let's solve for this numerically:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def chareq(x):\n",
    "    Kc, omega = x\n",
    "    s = 1j*omega\n",
    "    ce = 1 + Kc*loop\n",
    "    ce_of_s = ce(s)\n",
    "    return ce_of_s.real, ce_of_s.imag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.optimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([12.6       ,  0.89442719])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.optimize.fsolve(chareq, [13, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How can we determine stability without calculating the roots? See the [next notebook](SymPy%20Routh%20Array.ipynb)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}